function U_num = ex23_SolveSystem(bc_homog,N,id_bd,x_n,s_nears,p_s,dp_s,x_s,w_s)
%
% id_bd : indices for the boundary nodes with dirichlet BC's

sPts  = length(x_s);


%Stiffness matrix
K = sparse(zeros(N,N));
for k=1:sPts
  k_near = s_nears{k};
  dp_k   = dp_s{k};
  K(k_near,k_near) = K(k_near,k_near) + w_s(k)*dp_k*dp_k';
end


% right hand side computation
f  = zeros(N,1);
nB = length(id_bd);
% 0: homogeneous BC; 1: non-homogeneous BC
if bc_homog == 0
  for k=1:sPts
    k_near = s_nears{k};
    n_k    = length(k_near);
    p_k    = p_s{k};
    w_gk   = w_s(k);			
    val    = -2.0 * w_gk * x_s(k,2);
    % loop on the neighbors of the i-th gauss points
    for ia=1:n_k
      i    = k_near(ia);
      f(i) = f(i) + val*p_k(ia);
    end
  end
  for j=1:nB
    i    = id_bd(j);
    f(i) = 0;
  end
else
  % Non-Homogeneous Boundary Conditions
  for j=1:nB
    i   = id_bd(j);
    b   = x_n(i,2)*x_n(i,1)*(1-x_n(i,1));
    f(:)= f(:) - b*K(:,i);
  end
  for j=1:nB
    i    = id_bd(j);
    f(i) = x_n(i,2)*x_n(i,1)*(1-x_n(i,1));
  end
end

for j=1:nB
  i      = id_bd(j);
  K(i,:) = 0; 
  K(:,i) = 0;
  K(i,i) = 1;
end

% The solution is computed
U_num = K\f;